Complex SIR Model

Assumption

  • The recovered (R) compartment will no longer be susceptible

Parameters with default values:

  • Population \((N = 300)\)

  • Transmission/Infection Rate \((r = 0.02)\)

  • Percentage of healthy individuals with symptoms \((g = 0.2)\)

  • Initial number of infected individuals day1 \((d_1 = 5)\)

  • Initial number of cured individuals day1 \((c_1 = 0)\)

  • Percentage of individuals receiving Treatment A \((PTA = 1)\)

  • Percentage of individuals receiving Treatment B \((PTB = 0)\)

  • Treatment A effectiveness \((A_{eff} = 0.75)\)

  • Treatment B effectiveness \((B_{eff} = 0.25)\)

  • Treatment A cost \((A_{cost} = 20.00)\)

  • Treatment B cost \((B_{cost} = 15.00)\)

  • Number of days since first case \((day = 38)\)

Model upgrade with randomness - Additional Variables

  • Initial number of known cases (untreated) day1 \((u_1 = 0)\)

  • Percentage of individual being tested daily \((PS = 1)\)

  • Test sensitivity \((tp = 0.99)\)

  • Test specificity \((tn = 0.95)\)

  • Testing cost \((Test.cost = 1)\)

complex.model.rand = function(N=300, r=0.02, g=0.2, d1=5, c1=0, u1=0, PS=1, tp=0.99, tn=0.95, PTA=1, PTB=0, A.eff=0.75, B.eff=0.25, day = -1, A.cost = 20, B.cost = 15, Test.cost = 1){
  
  # Initial conditions
  curr.day = 1
  d = d1
  c = c1
  u = u1
  # Calculate initial statistics
  h = N - d - c
  hs = round(g * h)
  s = hs + d
  a = s - u
  ns = round(PS * a)
  ad = round(PS * (d - u))
  as = ns - ad
  tps = rbinom(1, as, 1-tn)
  tpd = rbinom(1, ad, tp)
  tpt = tps + tpd
  tpt1 = tpt
  NDAT = tpd + u
  NDTA = 0
  NDTB = 0
  DCA = 0
  DCB = 0
  NAT = 0
  NTA = 0
  NTB = 0
  SDC = 0
  New = min(N - c - d, round(r * hs * d))
  
  # Table storing S, I, R compartments over time
  df = data_frame(Day = curr.day, 
                  Healthy = h,
                  Symptopms = hs,
                  Diseased = d,
                  Tot.Symp = s,
                  Avail.screen = a,
                  Untreated = u,
                  Screen.Tot = ns,
                  Screen.Sym = as,
                  Screen.Dis = ad,
                  Test.Pos.Sym = tps,
                  Test.Pos.Dis = tpd,
                  Test.Pos.Total = tpt,
                  Avail.TRT.Dis = NDAT,
                  TRT.A.Dis = NDTA, 
                  TRT.B.Dis = NDTB, 
                  TRT.A.Cure = DCA, 
                  TRT.B.Cure = DCB,
                  Avail.TRT.Tot = NAT,
                  TRT.A.Tot = NTA,
                  TRT.B.Tot = NTB,
                  Sum.Dis.Cured = SDC,
                  Cured = c,
                  Catch.Dis = New)
  
  # Predict S, I, R compartments by day
  while(day == -1 & d > 0 | (day != -1 & curr.day < day)){
    # if number of days is specified, calculate up until that day
    # else calculate until there is no more diseased individual
    
    # Update current day and c, u, d
    curr.day = curr.day + 1
    c = c + SDC
    u = d - SDC
    d = d + New - SDC
    
    # Calculate other variables
    h = N - d - c
    hs = round(g * h)
    s = hs + d
    a = s - u
    ns = round(PS * a)
    ad = round(PS * (d - u))
    as = ns - ad
    tps = rbinom(1, as, 1-tn)
    tpd = rbinom(1, ad, tp)
    tpt = tps + tpd
    NDAT = tpd + u
    NDTA = round(PTA * NDAT)
    NDTB = min(NDAT - NDTA, round(PTB * NDAT))
    DCA = min(NDTA, round(NDTA * A.eff))
    DCB = min(NDTB, round(NDTB * B.eff))
    if (curr.day == 2){
      NAT = tpt + tpt1
    }
    else{
      NAT = tpt + u
    }
    NTA = round(PTA * NAT)
    NTB = NAT - NTA
    SDC = DCA + DCB
    New = min(N - c - d, round(r * hs * d))
    
    data = c(curr.day, h, hs, d, s, a, u,
             ns, as, ad, tps, tpd, tpt,
             NDAT, NDTA, NDTA, DCA, DCB, NAT, NTA, NTB,
             SDC, c, New)
    df = rbind(df, data)
  }
  df$Cost = df$TRT.A.Tot * A.cost + df$TRT.B.Tot * B.cost + df$Screen.Tot * Test.cost
  return(df)
}

model.plot = function(model){
  cost = sum(model$Cost)
  sick = sum(model$Diseased)
  peak = max(model$Diseased)
  ggplot(data  = model, aes(x = Day)) +
    geom_line(aes(y = Healthy, color = "Healthy")) +
    geom_line(aes(y = Diseased, color = "Diseased")) +
    geom_line(aes(y = Cured, color = "Cured")) +
    scale_color_discrete(name = "Compartment") +
    labs(y = "Population", title = paste("Complex SIR Model: infected = " , sick, ", peak = ", peak, ", cost = $", cost, sep = "")) +
    theme_minimal()
}
df = complex.model.rand()
df
model.plot(df)

Checking variability by running 100 replications of the model

line.plot = function(model){
  list(geom_point(data = model, aes(x = Day, y = Healthy, color = "Healthy"), size = 0.8),
       geom_point(data = model, aes(x = Day, y = Diseased, color = "Diseased"), size = 0.8),
       geom_point(data = model, aes(x = Day, y = Cured, color = "Cured"), size = 0.8),
       geom_line(data = model, aes(x = Day, y = Healthy, color = "Healthy"), size = 0.1),
       geom_line(data = model, aes(x = Day, y = Diseased, color = "Diseased"), size = 0.1),
       geom_line(data = model, aes(x = Day, y = Cured, color = "Cured"), size = 0.1))
}

model.rep.plot = function(reps){
  plot = ggplot() + labs(y = "Population", title = paste("Complex SIR Model with", reps, "reps")) + theme_minimal()
  for(i in 1:reps){
    df = complex.model.rand()
    plot = plot + line.plot(df)
  }
  plot = plot + scale_color_discrete(name = "Compartment")
  plot
  return(plot)
}
model.rep.plot(10)

model.rep.plot(100)

accumulate_by <- function(dat, var) {
  var <- lazyeval::f_eval(var, dat)
  lvls <- plotly:::getLevels(var)
  dats <- lapply(seq_along(lvls), function(x) {
    cbind(dat[var %in% lvls[seq(1, x)], ], frame = lvls[[x]])
  })
  dplyr::bind_rows(dats)
}
visual_data = select(df, c(Day, Healthy, Diseased, Cured)) %>% tidyr::gather(key = Compartment, value = Population, Healthy, Diseased, Cured) %>% accumulate_by(~Day)
plot_ly(data = visual_data, type = "scatter", mode = "lines",
        x = ~Day, y = ~Population, color = ~Compartment, frame = ~frame) %>% 
  layout(hovermode = 'compare') %>%
  animation_opts(frame = 100, transition = 0, redraw = FALSE, mode = "immediate") %>%
  animation_slider(hide = T) %>%
  animation_button(x = 1, xanchor = "right", y = -0.1, yanchor = "bottom", label = "View Plot")
`arrange_()` is deprecated as of dplyr 0.7.0.
Please use `arrange()` instead.
See vignette('programming') for more help
This warning is displayed once every 8 hours.
Call `lifecycle::last_warnings()` to see where this warning was generated.
LS0tCnRpdGxlOiAiQ29tcGxleCBTSVIgbW9kZWwgd2l0aCBWYXJpYWJpbGl0eSIKYXV0aG9yOiAiTGluaCwgQnJpdG5leSwgQm93ZW4iCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCmBgYHtyIHNldHVwLCBpbmNsdWRlPUZBTFNFfQprbml0cjo6b3B0c19jaHVuayRzZXQod2FybmluZyA9IFRSVUUsIG1lc3NhZ2UgPSBUUlVFKQpsaWJyYXJ5KG1vc2FpYykKbGlicmFyeShnZ3Bsb3QyKQpgYGAKCiMjIENvbXBsZXggU0lSIE1vZGVsCgojIyMgQXNzdW1wdGlvbgoKKiBUaGUgcmVjb3ZlcmVkIChSKSBjb21wYXJ0bWVudCB3aWxsIG5vIGxvbmdlciBiZSBzdXNjZXB0aWJsZQoKIyMjIFBhcmFtZXRlcnMgd2l0aCBkZWZhdWx0IHZhbHVlczoKCiogUG9wdWxhdGlvbiAkKE4gPSAzMDApJAoKKiBUcmFuc21pc3Npb24vSW5mZWN0aW9uIFJhdGUgJChyID0gMC4wMikkCgoqIFBlcmNlbnRhZ2Ugb2YgaGVhbHRoeSBpbmRpdmlkdWFscyB3aXRoIHN5bXB0b21zICQoZyA9IDAuMikkCgoqIEluaXRpYWwgbnVtYmVyIG9mIGluZmVjdGVkIGluZGl2aWR1YWxzIGRheTEgJChkXzEgPSA1KSQKCiogSW5pdGlhbCBudW1iZXIgb2YgY3VyZWQgaW5kaXZpZHVhbHMgZGF5MSAkKGNfMSA9IDApJAoKKiBQZXJjZW50YWdlIG9mIGluZGl2aWR1YWxzIHJlY2VpdmluZyBUcmVhdG1lbnQgQSAkKFBUQSA9IDEpJAoKKiBQZXJjZW50YWdlIG9mIGluZGl2aWR1YWxzIHJlY2VpdmluZyBUcmVhdG1lbnQgQiAkKFBUQiA9IDApJAoKKiBUcmVhdG1lbnQgQSBlZmZlY3RpdmVuZXNzICQoQV97ZWZmfSA9IDAuNzUpJAoKKiBUcmVhdG1lbnQgQiBlZmZlY3RpdmVuZXNzICQoQl97ZWZmfSA9IDAuMjUpJAoKKiBUcmVhdG1lbnQgQSBjb3N0ICQoQV97Y29zdH0gPSAyMC4wMCkkCgoqIFRyZWF0bWVudCBCIGNvc3QgJChCX3tjb3N0fSA9IDE1LjAwKSQKCiogTnVtYmVyIG9mIGRheXMgc2luY2UgZmlyc3QgY2FzZSAkKGRheSA9IDM4KSQKCioqTW9kZWwgdXBncmFkZSB3aXRoIHJhbmRvbW5lc3MgLSBBZGRpdGlvbmFsIFZhcmlhYmxlcyoqCgoqIEluaXRpYWwgbnVtYmVyIG9mIGtub3duIGNhc2VzICh1bnRyZWF0ZWQpIGRheTEgJCh1XzEgPSAwKSQKCiogUGVyY2VudGFnZSBvZiBpbmRpdmlkdWFsIGJlaW5nIHRlc3RlZCBkYWlseSAkKFBTID0gMSkkCgoqIFRlc3Qgc2Vuc2l0aXZpdHkgJCh0cCA9IDAuOTkpJAoKKiBUZXN0IHNwZWNpZmljaXR5ICQodG4gPSAwLjk1KSQKCiogVGVzdGluZyBjb3N0ICQoVGVzdC5jb3N0ID0gMSkkCgpgYGB7ciBtb2RlbC5yYW5kfQpjb21wbGV4Lm1vZGVsLnJhbmQgPSBmdW5jdGlvbihOPTMwMCwgcj0wLjAyLCBnPTAuMiwgZDE9NSwgYzE9MCwgdTE9MCwgUFM9MSwgdHA9MC45OSwgdG49MC45NSwgUFRBPTEsIFBUQj0wLCBBLmVmZj0wLjc1LCBCLmVmZj0wLjI1LCBkYXkgPSAtMSwgQS5jb3N0ID0gMjAsIEIuY29zdCA9IDE1LCBUZXN0LmNvc3QgPSAxKXsKICAKICAjIEluaXRpYWwgY29uZGl0aW9ucwogIGN1cnIuZGF5ID0gMQogIGQgPSBkMQogIGMgPSBjMQogIHUgPSB1MQogICMgQ2FsY3VsYXRlIGluaXRpYWwgc3RhdGlzdGljcwogIGggPSBOIC0gZCAtIGMKICBocyA9IHJvdW5kKGcgKiBoKQogIHMgPSBocyArIGQKICBhID0gcyAtIHUKICBucyA9IHJvdW5kKFBTICogYSkKICBhZCA9IHJvdW5kKFBTICogKGQgLSB1KSkKICBhcyA9IG5zIC0gYWQKICB0cHMgPSByYmlub20oMSwgYXMsIDEtdG4pCiAgdHBkID0gcmJpbm9tKDEsIGFkLCB0cCkKICB0cHQgPSB0cHMgKyB0cGQKICB0cHQxID0gdHB0CiAgTkRBVCA9IHRwZCArIHUKICBORFRBID0gMAogIE5EVEIgPSAwCiAgRENBID0gMAogIERDQiA9IDAKICBOQVQgPSAwCiAgTlRBID0gMAogIE5UQiA9IDAKICBTREMgPSAwCiAgTmV3ID0gbWluKE4gLSBjIC0gZCwgcm91bmQociAqIGhzICogZCkpCiAgCiAgIyBUYWJsZSBzdG9yaW5nIFMsIEksIFIgY29tcGFydG1lbnRzIG92ZXIgdGltZQogIGRmID0gZGF0YV9mcmFtZShEYXkgPSBjdXJyLmRheSwgCiAgICAgICAgICAgICAgICAgIEhlYWx0aHkgPSBoLAogICAgICAgICAgICAgICAgICBTeW1wdG9wbXMgPSBocywKICAgICAgICAgICAgICAgICAgRGlzZWFzZWQgPSBkLAogICAgICAgICAgICAgICAgICBUb3QuU3ltcCA9IHMsCiAgICAgICAgICAgICAgICAgIEF2YWlsLnNjcmVlbiA9IGEsCiAgICAgICAgICAgICAgICAgIFVudHJlYXRlZCA9IHUsCiAgICAgICAgICAgICAgICAgIFNjcmVlbi5Ub3QgPSBucywKICAgICAgICAgICAgICAgICAgU2NyZWVuLlN5bSA9IGFzLAogICAgICAgICAgICAgICAgICBTY3JlZW4uRGlzID0gYWQsCiAgICAgICAgICAgICAgICAgIFRlc3QuUG9zLlN5bSA9IHRwcywKICAgICAgICAgICAgICAgICAgVGVzdC5Qb3MuRGlzID0gdHBkLAogICAgICAgICAgICAgICAgICBUZXN0LlBvcy5Ub3RhbCA9IHRwdCwKICAgICAgICAgICAgICAgICAgQXZhaWwuVFJULkRpcyA9IE5EQVQsCiAgICAgICAgICAgICAgICAgIFRSVC5BLkRpcyA9IE5EVEEsIAogICAgICAgICAgICAgICAgICBUUlQuQi5EaXMgPSBORFRCLCAKICAgICAgICAgICAgICAgICAgVFJULkEuQ3VyZSA9IERDQSwgCiAgICAgICAgICAgICAgICAgIFRSVC5CLkN1cmUgPSBEQ0IsCiAgICAgICAgICAgICAgICAgIEF2YWlsLlRSVC5Ub3QgPSBOQVQsCiAgICAgICAgICAgICAgICAgIFRSVC5BLlRvdCA9IE5UQSwKICAgICAgICAgICAgICAgICAgVFJULkIuVG90ID0gTlRCLAogICAgICAgICAgICAgICAgICBTdW0uRGlzLkN1cmVkID0gU0RDLAogICAgICAgICAgICAgICAgICBDdXJlZCA9IGMsCiAgICAgICAgICAgICAgICAgIENhdGNoLkRpcyA9IE5ldykKICAKICAjIFByZWRpY3QgUywgSSwgUiBjb21wYXJ0bWVudHMgYnkgZGF5CiAgd2hpbGUoZGF5ID09IC0xICYgZCA+IDAgfCAoZGF5ICE9IC0xICYgY3Vyci5kYXkgPCBkYXkpKXsKICAgICMgaWYgbnVtYmVyIG9mIGRheXMgaXMgc3BlY2lmaWVkLCBjYWxjdWxhdGUgdXAgdW50aWwgdGhhdCBkYXkKICAgICMgZWxzZSBjYWxjdWxhdGUgdW50aWwgdGhlcmUgaXMgbm8gbW9yZSBkaXNlYXNlZCBpbmRpdmlkdWFsCiAgICAKICAgICMgVXBkYXRlIGN1cnJlbnQgZGF5IGFuZCBjLCB1LCBkCiAgICBjdXJyLmRheSA9IGN1cnIuZGF5ICsgMQogICAgYyA9IGMgKyBTREMKICAgIHUgPSBkIC0gU0RDCiAgICBkID0gZCArIE5ldyAtIFNEQwogICAgCiAgICAjIENhbGN1bGF0ZSBvdGhlciB2YXJpYWJsZXMKICAgIGggPSBOIC0gZCAtIGMKICAgIGhzID0gcm91bmQoZyAqIGgpCiAgICBzID0gaHMgKyBkCiAgICBhID0gcyAtIHUKICAgIG5zID0gcm91bmQoUFMgKiBhKQogICAgYWQgPSByb3VuZChQUyAqIChkIC0gdSkpCiAgICBhcyA9IG5zIC0gYWQKICAgIHRwcyA9IHJiaW5vbSgxLCBhcywgMS10bikKICAgIHRwZCA9IHJiaW5vbSgxLCBhZCwgdHApCiAgICB0cHQgPSB0cHMgKyB0cGQKICAgIE5EQVQgPSB0cGQgKyB1CiAgICBORFRBID0gcm91bmQoUFRBICogTkRBVCkKICAgIE5EVEIgPSBtaW4oTkRBVCAtIE5EVEEsIHJvdW5kKFBUQiAqIE5EQVQpKQogICAgRENBID0gbWluKE5EVEEsIHJvdW5kKE5EVEEgKiBBLmVmZikpCiAgICBEQ0IgPSBtaW4oTkRUQiwgcm91bmQoTkRUQiAqIEIuZWZmKSkKICAgIGlmIChjdXJyLmRheSA9PSAyKXsKICAgICAgTkFUID0gdHB0ICsgdHB0MQogICAgfQogICAgZWxzZXsKICAgICAgTkFUID0gdHB0ICsgdQogICAgfQogICAgTlRBID0gcm91bmQoUFRBICogTkFUKQogICAgTlRCID0gTkFUIC0gTlRBCiAgICBTREMgPSBEQ0EgKyBEQ0IKICAgIE5ldyA9IG1pbihOIC0gYyAtIGQsIHJvdW5kKHIgKiBocyAqIGQpKQogICAgCiAgICBkYXRhID0gYyhjdXJyLmRheSwgaCwgaHMsIGQsIHMsIGEsIHUsCiAgICAgICAgICAgICBucywgYXMsIGFkLCB0cHMsIHRwZCwgdHB0LAogICAgICAgICAgICAgTkRBVCwgTkRUQSwgTkRUQSwgRENBLCBEQ0IsIE5BVCwgTlRBLCBOVEIsCiAgICAgICAgICAgICBTREMsIGMsIE5ldykKICAgIGRmID0gcmJpbmQoZGYsIGRhdGEpCiAgfQogIGRmJENvc3QgPSBkZiRUUlQuQS5Ub3QgKiBBLmNvc3QgKyBkZiRUUlQuQi5Ub3QgKiBCLmNvc3QgKyBkZiRTY3JlZW4uVG90ICogVGVzdC5jb3N0CiAgcmV0dXJuKGRmKQp9Cgptb2RlbC5wbG90ID0gZnVuY3Rpb24obW9kZWwpewogIGNvc3QgPSBzdW0obW9kZWwkQ29zdCkKICBzaWNrID0gc3VtKG1vZGVsJERpc2Vhc2VkKQogIHBlYWsgPSBtYXgobW9kZWwkRGlzZWFzZWQpCiAgZ2dwbG90KGRhdGEgID0gbW9kZWwsIGFlcyh4ID0gRGF5KSkgKwogICAgZ2VvbV9saW5lKGFlcyh5ID0gSGVhbHRoeSwgY29sb3IgPSAiSGVhbHRoeSIpKSArCiAgICBnZW9tX2xpbmUoYWVzKHkgPSBEaXNlYXNlZCwgY29sb3IgPSAiRGlzZWFzZWQiKSkgKwogICAgZ2VvbV9saW5lKGFlcyh5ID0gQ3VyZWQsIGNvbG9yID0gIkN1cmVkIikpICsKICAgIHNjYWxlX2NvbG9yX2Rpc2NyZXRlKG5hbWUgPSAiQ29tcGFydG1lbnQiKSArCiAgICBsYWJzKHkgPSAiUG9wdWxhdGlvbiIsIHRpdGxlID0gcGFzdGUoIkNvbXBsZXggU0lSIE1vZGVsOiBpbmZlY3RlZCA9ICIgLCBzaWNrLCAiLCBwZWFrID0gIiwgcGVhaywgIiwgY29zdCA9ICQiLCBjb3N0LCBzZXAgPSAiIikpICsKICAgIHRoZW1lX21pbmltYWwoKQp9CgpgYGAKCmBgYHtyfQpkZiA9IGNvbXBsZXgubW9kZWwucmFuZCgpCmRmCm1vZGVsLnBsb3QoZGYpCmBgYAoKIyMjIENoZWNraW5nIHZhcmlhYmlsaXR5IGJ5IHJ1bm5pbmcgMTAwIHJlcGxpY2F0aW9ucyBvZiB0aGUgbW9kZWwKCmBgYHtyIHBsb3RfZnVuY30KbGluZS5wbG90ID0gZnVuY3Rpb24obW9kZWwpewogIGxpc3QoZ2VvbV9wb2ludChkYXRhID0gbW9kZWwsIGFlcyh4ID0gRGF5LCB5ID0gSGVhbHRoeSwgY29sb3IgPSAiSGVhbHRoeSIpLCBzaXplID0gMC44KSwKICAgICAgIGdlb21fcG9pbnQoZGF0YSA9IG1vZGVsLCBhZXMoeCA9IERheSwgeSA9IERpc2Vhc2VkLCBjb2xvciA9ICJEaXNlYXNlZCIpLCBzaXplID0gMC44KSwKICAgICAgIGdlb21fcG9pbnQoZGF0YSA9IG1vZGVsLCBhZXMoeCA9IERheSwgeSA9IEN1cmVkLCBjb2xvciA9ICJDdXJlZCIpLCBzaXplID0gMC44KSwKICAgICAgIGdlb21fbGluZShkYXRhID0gbW9kZWwsIGFlcyh4ID0gRGF5LCB5ID0gSGVhbHRoeSwgY29sb3IgPSAiSGVhbHRoeSIpLCBzaXplID0gMC4xKSwKICAgICAgIGdlb21fbGluZShkYXRhID0gbW9kZWwsIGFlcyh4ID0gRGF5LCB5ID0gRGlzZWFzZWQsIGNvbG9yID0gIkRpc2Vhc2VkIiksIHNpemUgPSAwLjEpLAogICAgICAgZ2VvbV9saW5lKGRhdGEgPSBtb2RlbCwgYWVzKHggPSBEYXksIHkgPSBDdXJlZCwgY29sb3IgPSAiQ3VyZWQiKSwgc2l6ZSA9IDAuMSkpCn0KCm1vZGVsLnJlcC5wbG90ID0gZnVuY3Rpb24ocmVwcyl7CiAgcGxvdCA9IGdncGxvdCgpICsgbGFicyh5ID0gIlBvcHVsYXRpb24iLCB0aXRsZSA9IHBhc3RlKCJDb21wbGV4IFNJUiBNb2RlbCB3aXRoIiwgcmVwcywgInJlcHMiKSkgKyB0aGVtZV9taW5pbWFsKCkKICBmb3IoaSBpbiAxOnJlcHMpewogICAgZGYgPSBjb21wbGV4Lm1vZGVsLnJhbmQoKQogICAgcGxvdCA9IHBsb3QgKyBsaW5lLnBsb3QoZGYpCiAgfQogIHBsb3QgPSBwbG90ICsgc2NhbGVfY29sb3JfZGlzY3JldGUobmFtZSA9ICJDb21wYXJ0bWVudCIpCiAgcGxvdAogIHJldHVybihwbG90KQp9CmBgYAoKYGBge3J9Cm1vZGVsLnJlcC5wbG90KDEwKQptb2RlbC5yZXAucGxvdCgxMDApCmBgYApgYGB7cn0KYWNjdW11bGF0ZV9ieSA8LSBmdW5jdGlvbihkYXQsIHZhcikgewogIHZhciA8LSBsYXp5ZXZhbDo6Zl9ldmFsKHZhciwgZGF0KQogIGx2bHMgPC0gcGxvdGx5Ojo6Z2V0TGV2ZWxzKHZhcikKICBkYXRzIDwtIGxhcHBseShzZXFfYWxvbmcobHZscyksIGZ1bmN0aW9uKHgpIHsKICAgIGNiaW5kKGRhdFt2YXIgJWluJSBsdmxzW3NlcSgxLCB4KV0sIF0sIGZyYW1lID0gbHZsc1tbeF1dKQogIH0pCiAgZHBseXI6OmJpbmRfcm93cyhkYXRzKQp9CnZpc3VhbF9kYXRhID0gc2VsZWN0KGRmLCBjKERheSwgSGVhbHRoeSwgRGlzZWFzZWQsIEN1cmVkKSkgJT4lIHRpZHlyOjpnYXRoZXIoa2V5ID0gQ29tcGFydG1lbnQsIHZhbHVlID0gUG9wdWxhdGlvbiwgSGVhbHRoeSwgRGlzZWFzZWQsIEN1cmVkKSAlPiUgYWNjdW11bGF0ZV9ieSh+RGF5KQpwbG90X2x5KGRhdGEgPSB2aXN1YWxfZGF0YSwgdHlwZSA9ICJzY2F0dGVyIiwgbW9kZSA9ICJsaW5lcyIsCiAgICAgICAgeCA9IH5EYXksIHkgPSB+UG9wdWxhdGlvbiwgY29sb3IgPSB+Q29tcGFydG1lbnQsIGZyYW1lID0gfmZyYW1lKSAlPiUgCiAgbGF5b3V0KGhvdmVybW9kZSA9ICdjb21wYXJlJykgJT4lCiAgYW5pbWF0aW9uX29wdHMoZnJhbWUgPSAxMDAsIHRyYW5zaXRpb24gPSAwLCByZWRyYXcgPSBGQUxTRSwgbW9kZSA9ICJpbW1lZGlhdGUiKSAlPiUKICBhbmltYXRpb25fc2xpZGVyKGhpZGUgPSBUKSAlPiUKICBhbmltYXRpb25fYnV0dG9uKHggPSAxLCB4YW5jaG9yID0gInJpZ2h0IiwgeSA9IC0wLjEsIHlhbmNob3IgPSAiYm90dG9tIiwgbGFiZWwgPSAiVmlldyBQbG90IikKYGBgCgo=